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All  Ozbek,  and  Bernard  C.  Levy,  University  of  California  at  Davis 


S  16.4 


Summary 

The  multidimensional  inverae  scattering  problem  for  an  acous¬ 
tic  medium  is  considered  within  the  homogeneous  background 
Born  approximation.  The  objective  is  to  reconstruct  simulta¬ 
neously  the  velocity  and  density  profiles  of  the  medium.  The 
medium  is  probed  by  wide- band  plane- wave  sources,  and  the  time 
traces  observed  at  the  receivers  are  appropriately  filtered  to  ob¬ 
tain  generalized  projections  of  the  velocity  and  density  scattering 
potentials,  which  are  related  to  the  velocity  and  density  varia¬ 
tions  in  the  medium.  The  generalized  projections  are  weighted 
integrals  of  the  scattering  potentials;  in  the  two-dimensional  ge¬ 
ometry  the  weighting  functions  are  concentrated  along  parabo¬ 
las.  The  reconstruction  problem  for  the  generalized  projections  is 
formulated  in  a  way  similar  to  the  problem  of  x-ray,  or  straight- 
line  tomography.  The  solution  is  expressed  as  a  backprojection 
operation  followed  by  a  two  dimensional  space-invariant  filter¬ 
ing  operation.  In  the  Fourier  domain,  the  resulting  image  is  a 
linear  combination  of  the  velocity  and  density  scattering  poten¬ 
tials,  where  the  coefficients  depend  on  the  angle  of  incidence  of 
the  probing  wave.  Therefore,  two  or  more  different  angles  of 
incidence  are  necessary  to  solve  for  the  velocity  and  density  scat¬ 
tering  potentiab  separately. 

The  technique  of  defining  a  backprojection  operator  and  relat¬ 
ing  it  to  the  unknown  medium  for  the  case  of  zero-offset  problems, 
where  projections  over  circles  arise,  was  introduced  by  Fawcett 
(1985).  With  a  similar  technique,  Ozbek  &  Levy  (1987)  solved 
the  velocity  inversion  problem  in  constant-density  acoustic  media 
under  plane-wave  illumination,  where  parabolic  projections  are 
the  data.  This  work  extends  this  work  to  the  joint  reconstruction 
of  velocity  eind  density.  Only  the  2D  case  is  presented  here,  for 
the  3D  case  and  more  detailed  development,  see  Ozbek  &  Levy 
(1988).  The  extension  of  these  results  to  the  elastic  problem, 
where  elliptic  and  hyperbolic,  as  well  as  parabolic  projections 
are  inverted,  will  be  presented  elsewhere. 

Introduction 

Consider  the  scattering  experiment  described  in  Fig.  1.  A 
2-D  acoustic  medium  is  probed  by  a  wide-band  plane  wave  and 
the  scattered  field  is  observed  along  a  straight-line  receiver  array. 
The  pressure  field  P{x,u)  at  position  £=  (*,  y)  satisfies 

where  c(z)  is  the  propagation  velocity,  and  p(z)  is  the  density 
of  the  medium  at  point  Within  the  Born  approximation,  the 
scattered  field  Fi(if,b;)  at  receiver  location  (can  then  be  written 
as 

-  I^,(£')]i’.{£',a»)G.(i.s'.<v) 

+tf.(£')Vj.P.(x', w)  •  V^.G.(^,  x',  w)}  ,  (2) 


difference  between  the  total  field  P  and  the  incident  field  P^.  Go 
is  the  Green’s  function  associated  with  a  point  source  in  a  homo¬ 
geneous  medium.  Ue(s)  =  [cj/c*(£)]-  1  and  l/o(s)  =  ln[/>(i)//>„] 
can  be  termed  the  "velocity  scattering  potential"  and  the  "den¬ 
sity  scattering  potential",  respectively.  Co  is  the  propagation  ve¬ 
locity,  and  Po  is  the  density  of  the  background  medium.  We 
assume  that  c(x)  and  p(^  do  not  deviate  significantly  from  the 
background  values  of  Co  and  Po]  consequently  the  Uds)  and  ff,,(x) 
values  are  small  with  respect  to  1.  We  also  assume  that  Ife(x) 
and  have  the  bounded  support  V,  which  is  disconnected 

from  the  receiver  array. 

From  (2),  it  can  be  deduced  that  the  scattering  pattern  due  to 
ffc(£)  that  of  a  monopole,  whereas  the  scattering  pattern  due 
to  Uo{s)iB  that  of  the  sum  of  a  monopole  and  a  dipole.  Therefore, 
the  scattering  due  to  density  perturbations  is  most  prominent  for 
reflected  waves,  and  the  least  prominent  for  reflected  ones. 

The  incident  wave  is  given  as  Po(£',  w)  =  where  6  = 

(cos  9,  sin  9)  is  the  unit  vector  which  indicates  the  angle  of  in¬ 
cidence  of  the  plane-wave  source.  For  the  2-D  geometry,  the 
Green's  function  is  given  as  Go(x,i',w)  =  «H^**{lk|i  - 1'|)/4, 
where  Hg*^(')  indicates  the  Hankel  function  of  order  zero  and 
type  one. 

We  now  now  filter  the  observed  scattered  field: 

=  (3) 

From  (2),  the  inverse  Fourier  transform  of  §((,  Jfe,)  with  respect 
to  kr  can  be  written  as 


This  equation  expresses  ff(f,r)  as  a  weighted  integral  of  scatter¬ 
ing  potentiab  tfe(s)  and  Uf(s),  where  the  weighting  function  b 
nonzero  in  a  region  with  parabolic  support.  In  the  following,  it 
wUl  be  assumed  that  the  projections  j(f,r)  constitute  the  data 
that  b  given  by  the  scattering  experiment.  Therefore,  the  inverse 
scattering  problem  can  be  formulated  as  follows:  given  the  gen- 
eralbed  projections  {j((,r)  :  — oo  <(<oo,  0<r<  oo},  we 
want  to  reconstruct  the  scattering  potentbb  Cfe(x)  and  ff^(x). 

It  b  interesting  to  note  that  the  parabolical  projections  y  ((,  r) 
can  be  obtained  in  the  time  domain  abo: 

9{Cr)  -  -2rcof^  dr  f  dsP,((,s) 

CO  J —oo 

-b4x*Ce/  dr  f  dsP,(|,s).  (5) 

J—oo  J—oo 

The  Backprojection  Operation 


where  k  =  nr/cg  b  the  wavenumber.  The  scattered  field  P,  b  the 


The  first  step  of  our  inversion  procedure  b  to  perform  a  back- 
projection  operation  on  the  projections  y(f,r).  We  define  it  as 


2 


Velocity  and  density  inversion 


J-co  J-ao 


l(r  —  £•  I  —  |i—  Ci) 


J-oo  J-OO  y(r_£.l)2  _  |j_£|2 

At  a  given  point  £,  this  operation  sums  the  contributions  of  the 
projections  ff(f,r)  which  correspond  to  scattering  equation  (4). 
By  performing  this  backprojection  operation  for  every  point  in 
the  plane,  this  gives  an  image,  !7b(i). 

Our  first  objective  is  to  relate  s(f,r)  and  Ub{^  in  the  fre¬ 
quency  domain.  It  can  be  shown  that  the  Fourier  transform  of 
Ub{s)  is  given  by  (Ozbek  &  Levy  (1987)) 


Uam  =  fc  =  it,  = 

k-e  \  ^  2k-9  2k-9 


■where  g{k(,k,)  is  the  2-D  Fourier  transform  of  j(f,r),  fc  =  {kx,k„), 
k=  |k|,  A  =  (fc®  —  fc*,  2kxky),  ^  =  (cos(fl  +  0),  sin(9  -f-  ^)),  and 
^  =  (sin(5  -t- 1^),  —  cos(fl  -f-  ^)). 

Separate  Reconstruction  of  and 

In  this  section,  we  first  derive  a  frequency  domain  relation¬ 
ship  between  the  projections  9(f,>‘)  and  the  scattering  potentials 
I7c(£)  and  tfp(s),  thus  obtaining  a  "Projection  Slice  Theorem” 
eissociated  with  the  problem.  Inverting  this  relationship  provides 
both  a  frequency  domain  relationship  between  Uc{k),  and 

and  a  method  for  the  separate  reconstruction  of  I/c(x) 

and 

From  (2)  and  (3),  we  obtain  (Ozbek  As  Levy  (1988)) 

Hk(,kr)  =  |()7e  -  Up){k  =  kri+k(l^  + 

_  S  Up{k  =  kr9  +  k(±^  +  L  (8) 

for  |fcj|  <  |k,|,  where  Uc{k)  and  &,,(fc)  are  the  2-D  Fourier  trans¬ 
forms  of  i7a(i)  and  Up{^  respectively,  B  =  'ifsgn(kr)'\/fc?  - 


a  /  -1-1  ifi-^  —  p>0 
1  if  X  •  ^  —  p  <  0 


for  all  I  €  T, 
for  all  X  e  "V. 


For  |kj|  >  |k,|,  g(k(,kr)  is  related  to  the  part  of  the  observed 
scattered  field  that  corresponds  to  evanescent  waves  (Ozbek  ic 
Levy  (1987)),  and  we  do  not  make  use  of  this  portion  of  g{k(,  k,) 
in  our  inversion  scheme  .  The  inverse  formula  of  (8)  is 


Unm  = 


=  !2^e-pA-i/2*-£S 


2k-9’  ’’  2k- 9 


for  6  C,  where  the  cone  C  is  defined  below.  Uii{k)  denotes 
the  2-D  Fourier  transform  of  the  reconstructed  potential  Ub{^ 
which  is  obtained  by  applying  the  constant  density  reconstruc¬ 
tion  procedure  to  the  projections  s(f,r)  obtained  from  a  variable 
density  and  velocity  medium. 

Equation  (8)  represents  the  "Projection  Slice  Theorem”  asso¬ 
ciated  with  the  variable  density  inverse  acoustic  problem  relating 
the  1-D  Fourier  transform  of  s(f ,  k,)  with  respect  to  £  to  a  semi¬ 
circular  slice  of  the  2-D  Fourier  transform  of  CfR(i).  For  a  fixed 


^ri  9{k(>  kr)  gives  Ur{!^  along  a  semicircle  of  radius  |fcr|  centered 
at  fcr£.  By  letting  k,  vary,  these  semicircles  span  a  cone  C,  which 
is  defined  as 

O  =  {fc  ;  |fc  •  £1  >  \/2/2}  (10) 

for  7  =  -1-1,  where  k  =  {k,k)  and  £=  (cos(9-)-^)/2,  sin(fl-|-<^)/2). 
The  angular  range  of  this  cone  is  90°.  For  7  =  -1,  C  is  the 
complement  C  of  the  above  cone. 

Combining  (7)  and  (9)  gives 

UnUi)=^^UB(ki  =  U,{k)-2cos^(Up(ki,  keC.  (11) 

where  f  is  the  angle  between  the  vectors  k  and  £,  with  k  = 
(^k).  This  relation  shows  that  Ur{!^  can  be  computed  by  two- 
dimensional  filtering  of  the  backprojected  image  Uaik),  but  then 
since  Dr (k)  is  a  linear  combination  of  Ux{k)  and  Up{!^  we  cannot 
reconstruct  these  two  potentials  separately  from  a  single  exper¬ 
iment.  Therefore,  to  reconstruct  De(k)  and  D,,(fc)  separately,  in 
principle  we  need  two  experiments  with  plane  waves  incident  at 
angles  £1  and  £2;  then  we  can  solve  the  sy stein 

■  1  -2(fc-£i)^  1  f  1  =  f  fl2l 

1  -2(fc-£2)^J  [  D,(fc)  DR2(a  j  ’  ^  ^ 

Mith.L) 

which  requires  inverting  the  matrix  M{^9^,^). 

For  the  numerical  stability  and  robustness  of  the  matrix  in¬ 
version  procedure,  the  matrix  Af’(^£l,£2)  must  be  as  nonsingu¬ 
lar  as  possible.  The  most  appropriate  measure  of  the  singularity 
of  a  matrix  is  the  smallest  singular  value  of  the  matrix.  Inver¬ 
sion  of  M  would  be  most  robust  when  the  smallest  singular  value 
0'min(  Af )  is  maximized.  This  takes  place  for  values  of  fc,  £j ,  and  £2 
such  that  £j  •  £2  =  0  and  fc  =  ±£i  or  fc  =  ±£2.  Therefore  the  two 
probing  waves  are  incident  at  angles  perpendicular  to  each  other. 
Under  this  condition,  let  us  consider  the  frequency  domain  cov¬ 
erage  we  would  have  for  finite  bandwidth  data,  assuming  that  we 
have  receiver  coverage  surrounding  the  medium.  Neglecting  the 
low  frequency  cutofi'band,  the  frequency  domain  coverage  due  to 
a  single  probing  wave  has  a  "figure-of-eight”  shape  aligned  with 
the  direction  of  the  probing  wave  (Ozbek  ic  Levy  (1987)).  When 
two  probing  waves  are  used,  Ue{^  and  D,(k)  can  be  solved  only 
in  regions  where  where  there  is  double  coverage,  as  indicated  by 
shaded  areas  in  Fig.  2.  However,  if  we  consider  the  superimposed 
"radiation  pattern"  of  cr,nin(Af)  drawn  in  Fig.  2  also,  we  see  that 
M  is  most  singular  for  those  values  of  fc  where  we  have  double 
coverage. 

For  general  values  of  £j  and  £j  the  situation  is  similar.  In  gen¬ 
eral,  M  is  singular  for  values  of  fc  which  satisfy  [fc  •  £i|  =  |fc  •  £2|- 
Therefore,  for  £1  =  ±£2,  Af  is  singular  for  all  fc;  otherwise,  it  is 
singular  for  fc  =  ±(£j  +  £2!  or  i  =  ±{£i  -  £2)/|£l  “  ^1- 

These  are  the  directions  which  in  fact  bisect  the  regions  where 
there  is  double  coverage. 

In  practice,  then,  it  would  be  appropriate  to  use  more  than 

two  angles,  say  anglra  £l,  £2 . 9jf,  and  for  each  fc,  solve  the 

resulting  system 


1  -2(fc  ■£,..)’ 

1  -2(i-£.y)* 

1 

M{k) 


^Hu.(fc) 

UrhM 

ie(fc) 


Velocity  and  density  inversion 


3 


] 


by  the  least  squares  method,  where  {ty ,  «y, . . .  ,ttp} 

C  {1,2, . . . ,  JV}  is  the  set  of  indices  corresponding  to  the  angles 
of  incidence  for  which  the  probing  wave  provides  coverage  at  k. 
This  gives  the  solution 


UM 


(14) 


Numerical  Example 


The  theory  presented  in  this  paper  was  tested  for  the  two- 
dimensional  case,  using  computer-generated  synthetic  data.  Figs. 
3a  and  3b  show  the  velocity  and  density  scattering  potential  mod¬ 
els,  Uc{^  and  respectively.  The  scattering  potentials  cor¬ 

respond  to  velocity  and  density  anomalies  which  are  constant  in 
square-shaped  areas  of  dimensions  35  m  X  35  m.  The  background 
medium  was  homogeneous  with  velocity  5000  m/s.  The  source 
wavelet  wm  lowpass  with  a  cutoff  frequency  of  425  Hz,  so  that 
the  object  sizes  are  three  times  the  shortest  wavelength  in  the 
source  signal.  The  regions  of  anomaly  are  separated  by  a  dis¬ 
tance  six  times  the  shortest  wavelength.  The  synthetic  scattered 
waves  were  obtained  by  using  the  forward  scattering  equation 
under  Born  approximation;  however,  since  the  object  sizes  were 
not  too  large  with  respect  to  the  shortest  wavelength,  we  do  not 
deem  this  approximation  to  be  critical  for  this  example,  for  small 
scattering  potential  magnitudes.  The  entire  image  area  was  500 
m  X  500  m,  the  grid  size  was  5  m  x  5  m,  and  receivers  were 
located  on  all  sides  around  the  medium,  100  on  each  side. 

As  indicated  above,  for  numerical  stability  in  the  individ¬ 
ual  reconstruction  of  velocity  and  density  inhomogeneities,  more 
than  two  sources  are  needed.  In  this  experiment,  we  have  used 
eight  angles  of  incidence,  at  22.5”  intervals.  The  inversion  was 
performed  over  the  regions  in  the  k  domain  where  coverage  was 
provided  by  at  least  five  probing  waves;  i.e.,  using  the  notation  of 
eq.  (14),  iV  =  8,  P  >  5,  and  rank(M)  >  4  for  all  inversion  points 
k.  This  corresponds  to  carrying  out  the  inversion  over  a  circu¬ 
lar  lowpass  region  with  a  radius  of  about  55%  of  the  maximum 
frequency  coverage  provided  by  a  single  source. 

For  comparison,  we  first  present  the  images  Ub{^  and  Ur(x). 
To  obtain  these  images,  in  the  frequency  domain,  values  obtained 
due  to  different  sources  providing  multiple  coverage  were  simply 
averaged  point  by  point.  Fig.  4  shows  the  backprojected  image 
Ub{x).  Ub{3^  can  be  interpreted  as  a  "migrated”  image  of  the  ve¬ 
locity  field  for  a  constant  density  medium.  Fig.  5  depicts  Ur{x), 
which  is  the  image  obtained  by  applying  the  constant  density 
reconstruction  procedure  to  the  data  obtained  from  a  variable 
density  and  velocity  medium. 

Some  observations  can  be  made  regarding  these  images.  Both 
images  display  the  locations  of  the  scatterers;  however  the  "in¬ 
version”  image  looks  much  better  focused  than  the  "migration” 
image.  This  effect  in  general  has  been  noted  by  other  researchers 
also  (Esmersoy  &  Miller  (1987)).  In  addition,  the  values  of  i7B(i) 
are  orders  of  magnitude  different  than  the  original  scattering  po¬ 
tential  levels.  On  the  other  hand,  f7jj(z)  looks  like  !7c(®)  -  i7,(i), 
and  actual  constructed  values  confirm  this.  To  see  how  this  comes 
about,  consider  that  our  averaging  scheme  in  the  frequency  do¬ 
main  can  be  ideally  written  as 

Uriki  =  -  r  =  ^c(fc)  -  U.ik).  (15) 

JT  ^0 

using  eqn.  (9),  where  Ur  corresponds  to  the  inversion  result  for 
one  source,  and  Ur  corresponds  to  the  result  after  averaging. 


Therefore,  what  we  have  actually  obtained  this  way  is  the  scatter¬ 
ing  potential  associated  with  the  compressibility  of  the  medium. 

Figs.  6a  and  6b  show  the  separate  reconstructions  of  i/c(x) 
and  f/,,(x),  respectively.  The  numerical  values  obtained  are  within 
20%  of  the  model  values.  The  sources  of  error  are  the  bilinear 
interpolation  used,  the  finite  size  of  receiver  arrays  on  each  side, 
lowpass  nature  of  the  frequency  domain  coverage,  the  source  de- 
convolution  process,  and  the  lack  of  zero  frequency  information. 

In  fact,  the  the  scattered  field  P,(f,  w)  has  zero  amplitude  for 
zero  frequency,  as  one  can  observe  from  the  Lippmann-Schwinger 
equation.  In  addition,  the  coefficient  of  Uf(k)  in  eqn.  (9)  is  not 
analytic  around  k_  =  Q.-  Therefore,  the  DC  level  cannot  be  re¬ 
constructed  with  the  derived  inversion  formulas.  In  our  imple¬ 
mentation,  we  have  instead  estimated  it  from  the  closest  eight 
values,  arbitrarily  assigning  a  weight  of  1/6  to  the  closest  four 
samples,  and  1/12  to  the  next  closest  samples  which  are  diago¬ 
nally  located. 


Conclusions 

We  considered  the  problem  of  the  separate  reconstruction 
of  the  velocity  and  density  inhomogeneities  for  a  multidimen¬ 
sional  acoustic  medium  probed  by  wide-band  plane  waves.  The 
problem  was  posed  as  a  generalized  tomographic  problem,  where 
weighted  integrals  of  the  velocity  scattering  potential  Uc{ss)  and 
the  density  scattering  potential  Up{x)  are  considered  as  data.  A 
backprojection  operator  Ub(x)  was  defined,  which  was  related  to 
the  generalized  projections  in  the  Fourier  transform  domain.  It 
was  shown  that,  by  applying  a  time-invariant  filter  to  17jB(®),  we 
can  obtain  an  image,  Ur{x),  which  in  the  Fourier  domain  is  a 
linear  combination  of  the  velocity  and  density  scattering  poten¬ 
tials,  and  where  the  coefficients  depend  on  the  angle  of  incidence 
of  the  probing  wave.  Therefore,  for  numerical  stability,  several 
angles  of  incidence  were  used  to  solve  for  the  velocity  and  density 
scattering  potentials  separately. 

Acknowledgements 

This  work  was  supported  by  the  Air  Force  Office  of  Scientific 
Research  under  Grant  No.  AFOSR-85-0227,  by  the  National  Sci¬ 
ence  Foundation  under  Grant  No.  ECS-83-12921,  and  by  the  Ver¬ 
tical  Seismic  Profiling  Consortium  at  the  MIT  Earth  Resources 
Laboratory. 


References 

C.  Esmersoy  and  D.  Miller,  "Stacking  Versus  Back  Propaga¬ 
tion  in  Seismic  Imaging:  Duality  for  Multidimensional  Lin¬ 
earized  Inversion,”  presented  at  the  57th  Ann.  Int.  Mtg. 
and  Expos.,  Soc.  Explor.  Geophys.,  New  Orleans,  Expanded 
Abstracts,  pp.  744-746,  Oct.  1987. 

J.  A.  Fawcett,  "Inversion  of  N-Dimensional  Spherical  Averages,” 
SIAM  J.  Appl.  Math.,  Vol.  45,  No.  2,  pp.  336-341,  April 
1985. 

A.  Ozbek  and  B.  C.  Levy,  "Inversion  of  Parabolic  and  Paraboloid¬ 
al  Projections,”  submitted  to  IEEE  Trans.  Acoaat.  Speech 
Signal  Processing,  April  1987. 

A.  Ozbek  and  B.  C.  Levy,  "Simultaneous  Inversion  of  Velocity 
and  Density  Profiles  for  Multidimensional  Acoustic  Media,” 
submitted  to  J.  Aeoust.  Soc.  Am.,  May  1988. 


1083 


RegioBS  of 
double  coverage 


Fig.  2.  Frequency  coverage  of  UR(/f)  and  the  “radiation 
pattern”  of  (M;k:0i,  ^z)- 


